function [y] = spec_sub(x,Fs,l, event, autoflag, sil_ratio)
%  __  __                            _   _       
% |  \/  |                          | | | |      
% | \  / | __ _ _ __ __ _ _   _  ___| |_| |_ ___ 
% | |\/| |/ _` | '__/ _` | | | |/ _ \ __| __/ _ \
% | |  | | (_| | | | (_| | |_| |  __/ |_| ||  __/
% |_|  |_|\__,_|_|  \__, |\__,_|\___|\__|\__\___|
%                      | |                       
%                      |_|                       
%  _    _       _                    _ _         
% | |  | |     (_)                  (_) |        
% | |  | |_ __  ___   _____ _ __ ___ _| |_ _   _ 
% | |  | | '_ \| \ \ / / _ \ '__/ __| | __| | | |
% | |__| | | | | |\ V /  __/ |  \__ \ | |_| |_| |
%  \____/|_| |_|_| \_/ \___|_|  |___/_|\__|\__, |
%                                           __/ |
%                                          |___/ 
% Last Modified by Dexter Delfrate for 
% Marquette University's Speech and Signal Processing Lab.
%
% This file is part of the Marquette University filters add-on for XBAT.
%
% To run this program, open xbat.m and select 'Spectral Subtraction' from filters/Marquette.
% 
% Spectral subtraction is the oldest and simplest technique for signal
%   enhancement. In this method, the magnitude spectrum of the noise, estimated 
%   from a silence region, is literally subtracted from the magnitude spectrum 
%   of each frame, then recombined using the noisy phase spectrum to get the 
%   enhanced time domain signal.
% 
% Inputs:
%       x= noisy signal
%       Fs = sampling rate of the signal
%       l = framesize in points (default 256)
%       event = selection event defining the start and end times of a
%           silence region
%       autoflag = flag for auto silence-detect, 1=use (default 1)
%       sil_ratio = if autoflag=1, silence ratio (default 0.05)
%                    if autoflag=0, noise waveform to use (No default)
%              
% Outputs:
%       y = enhanced signal
%
%   See also
% 
%   Copyright 2009 Marquette University, Speech and Signal Processing
%   Laboratory.


framesize = l;
duration = event.duration;
times = event.time;
overlap = l/2;
frame_start = ceil(times(1)*Fs/framesize +1);
% frame_end =  floor(times(2)*Fs/128 +1);
x=x(:);
n=length(x);
x(n+1:n+4*l)=zeros(4*l,1);

num_frames=floor(n/(framesize-overlap)+1); %possibly should be framesize*overlap instead of minus?

%Noise spectrum estimation from the silence frames of data.
%disp('Noise Estimation...');
%noise1=x(1:l); 
%noise2=x(l+1:2*l);
%noise3=x(2*l+1:3*l);
%noise(1,:)=abs(fft(noise1'));
%noise(2,:)=abs(fft(noise2'));
%noise(3,:)=abs(fft(noise3'));
%ave_noise=mean(noise);
if autoflag==0
    nf = floor(duration*Fs/framesize +1); %number of noise frames
    for i=1:nf%num_win_noise
        An(:,i)=x(frame_start+(i-1)*l/2:frame_start + l-1 + (i-1)*l/2);
        %    An(:,i)=x(1+(i-1)*l/2:l + (i-1)*l/2);
    end
    N = abs(fft(An.*repmat(hanning(l),1,nf)));%num_win_noise)));
    clear An
    ave_noise=mean(N');
else
    silfr = findsil(x, framesize-overlap, sil_ratio);
    nf = length(silfr/2);
    for i=1:nf
        An(:,i)=x((1+(silfr(i)-1)*(framesize-overlap)):((silfr(i)+1)*(framesize-overlap)));
    end
    N = abs(fft(An.*repmat(hanning(l),1,nf)));%num_win_noise)));
    clear An
    ave_noise=mean(N');
end



% Computation of FFT(magnitude and phase) for each window
%disp('FFT Computation...');
for i=1:num_frames
   frame=x((i-1)*(l-overlap)+1:(i-1)*(l-overlap)+l);
   frame=frame.*hanning(l); %Triangular windowing
   frame_fft=fft(frame);
   ang_frame(:,i)=angle(frame_fft);
   mag_frame(:,i)=abs(frame_fft);
end

% Magnitude Averaging over 3 windows
%disp('Magnitude Averaging...');
for i=2:num_frames-1
   mag_frame(:,i)=(mag_frame(:,i-1)+mag_frame(:,i)+mag_frame(:,i+1))/3;
end

% Spectral Subtraction with Half-Wave rectification
%disp('Spectral Subtraction...');
for i=1:num_frames
   enhanced_mag=mag_frame(:,i)-ave_noise';
   for j=1:l
       if enhanced_mag(j) < 0
          enhanced_mag(j)=0;
       end
   end
% Phase information from noisy speech is added to enhanced speech
   fft_y=enhanced_mag.*cos(ang_frame(:,i))+j*enhanced_mag.*sin(ang_frame(:,i));
   time_y(:,i)=real(ifft(fft_y));
end

% Reconstruction(Hanning)
%disp('Reconstruction...');
y(1:l-overlap)=time_y(1:l-overlap,1);
for i=1:num_frames-1
   y(i*(l-overlap)+1:i*(l-overlap)+l-overlap)=time_y(l-overlap+1:l,i)+time_y(1:l-overlap,i+1);
end
y=y(1:n);
%subplot(211),plot(x(1:n));
%title('Noisy Speech Signal');
%ylabel('Amplitude');
%xlabel('Samples');
%ax=axis;
y=real(max(x)/max(real(y))*y);
y=y';
%subplot(212),plot(y);
%plot(y);
%title('Enhanced Speech Signal');
%ylabel('Amplitude');
%xlabel('Samples');
%axis(ax);

